PROGRAM MAPS
USE POLYMORPHIC_COMPLEXTAYLOR
TYPE(DAMAP) Id,Map,Map_pb
REAL(DP) L1,KICK1,L2,KICK2,LD 
REAL(DP) :: PREC = 1.D-10
type(vecfield) F
type(pbfield) H
integer mf,me
mf=20
open(unit=mf,file='map_of_fodo.txt')
CALL INIT(NO1=3,ND1=2,NP1=0,NDPT1 =0)     !   <------------------ init for maps in ND1 degrees of freedom


call alloc(Id,Map,Map_pb)
call alloc(F); call alloc(h);

KICK1=1.d0;   L1=0.5D0
KICK2=-0.9d0; L2=0.5D0
LD=1.5D0

WRITE(MF,*) " The beam line Q1-LD-Q2-LD is created. Remember Lie operators act in the reverse direction "

id=1
MAP=ID

call VECFIELD_QUAD(LD,0.D0,F); MAP=TEXP(F,MAP);
call VECFIELD_QUAD(L2,KICK2,F); MAP=TEXP(F,MAP);
call VECFIELD_QUAD(LD,0.D0,F); MAP=TEXP(F,MAP);
call VECFIELD_QUAD(L1,KICK1,F); MAP=TEXP(F,MAP);

CALL PRINT(MAP,MF)

WRITE(MF,*) " SAME MAP CREATED USING THE POISSON BRACKET OPERATOR "

id=1
Map_pb=ID

call PBFIELD_QUAD(LD,0.D0,H); Map_pb=TEXP(H,Map_pb);
call PBFIELD_QUAD(L2,KICK2,H); Map_pb=TEXP(H,Map_pb);
call PBFIELD_QUAD(LD,0.D0,H); Map_pb=TEXP(H,Map_pb);
call PBFIELD_QUAD(L1,KICK1,H); Map_pb=TEXP(H,Map_pb);

CALL PRINT(Map_pb,MF)

WRITE(MF,*) " TESTING AND COMPARING BOTH MAPS: ANSWER SHOULD BE IDENTITY TO ORDER NO "

id=map*Map_pb**(-1)
!ID=ID.CUT.C_%NO
CALL PRINT(ID,MFREC)

call kill(Id,Map,Map_pb)
call kill(F); call kill(h);
close(mf);
END PROGRAM MAPS

subroutine VECFIELD_QUAD(L,K,F)
USE POLYMORPHIC_COMPLEXTAYLOR
implicit none
TYPE(VECFIELD) F
REAL(DP) L,K

! DRIFT PART
F%V(1)=L*(1.D0.MONO.'01')  /SQRT(1.D0-(1.D0.MONO.'02')-(1.D0.MONO.'0002'))
F%V(3)=L*(1.D0.MONO.'0001')/SQRT(1.D0-(1.D0.MONO.'02')-(1.D0.MONO.'0002'))

! QUADRUPOLE PART 

F%V(2)=-L*K*(1.D0.MONO.'1')
F%V(4)= L*K*(1.D0.MONO.'001')


end subroutine VECFIELD_QUAD

subroutine PBFIELD_QUAD(L,K,H)
USE POLYMORPHIC_COMPLEXTAYLOR
implicit none
TYPE(PBFIELD) H
REAL(DP) L,K

! DRIFT PART
H=L*SQRT(1.D0-(1.D0.MONO.'02')-(1.D0.MONO.'0002'))

! QUADRUPOLE PART 

H%H=H%H-(L*K/2.D0)*((1.D0.MONO.'2')-(1.D0.MONO.'002'))

end subroutine PBFIELD_QUAD

